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Abstract 

By extending the nonequilibrium potential refinement algorithm(NEPR) and lattice switch(LW) 
method to the semigrand ensemble, the semigrand potentials of the fee and hep structure of 
polydisperse hard sphere crystals are calculated with the bias sampling scheme. The result shows 
that the fee structure is more stable than the hep structure for polydisperse hard sphere crystals 
below the terminal polydispersity. 
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I. INTRODUCTION 

A set of hard spheres under thermal agitation constitute a simple yet non trivial model of 
condensed matter and especially represents an idealization of a very important class of real 
colloid dispersions. The model has been extensively studied in the past decades. One of the 
important feature of the model is that the system undergoes a purely entropy-driven first- 
order phase transition from the fluid phase to a crystal phase at sufficiently high densityi^. 
Simple estimations of the free energy of different crystal structures reveal that the possible 
structure could be face-centered cubic (fee) or hexagonal close packed (hep). However, 
due to the very similarity in local environments of the two structures, the difference of free 
energy between them is extremely small and very hard to determine. The determination 
of the relative stability between the two structures from theoretical calculations has a long 
history 4 . A clear consensus was reached in the last decade that fee is the more stable 
phase^^, however, there are still different views on the problem. The recent results of 
Pronk and Frenkel^ indicate that a moderate deformation of a hard-sphere crystal may 
make the hep phase more stable than the fee phase. Kwak and Kofke^ investigated the 
effect of monovacancies on the relative stability of fee and hep hard-sphere crystals. 

The particle size of most artificial colloidal systems are polydisperse, the polydispersity 
is defined as the ratio of the standard deviation and the mean of the diameter distribution 
of particles, which is an intrinsic property of a colloidal system. The polydispersity may 
significantly affect the thermodynamic and dynamic properties of a hard sphere system, e.g., 
there exists a terminal polydispersity above which no cystallization can occur— iii^, the os- 
motic pressure of a polydisperse hard sphere crystal is higher than the one of a monodisperse 
system with the same volume fraction 1 ^ 4 -, and there are local fractionations of particle sizes 
which has a strong retarding effect on nucleation 1 ^^. However, the influence of size polydis- 
persity on the relative stability of fee and hep hard-sphere crystals has not been addressed. 
In this paper we will compute the free energy difference between polydisperse fee and hep 
hard sphere crystals by Monte Carlo (MC) simulations, which suggests that the fee phase is 
still more stable than hep phase below the terminal polydispersity. 

The simple simulating method in canonical ensemble is not suitable for a polydisperse 
system with continuous distribution of particle sizes. The reason is very simple, the simu- 
lation system is often too small to realize a given particle size distribution to the designed 



accuracy. Thus the grand canonical ensemble or semigrandr^ ensemble must be used. In 
these ensembles the number of particles of each size (characterized by the particle diameter 
<j), N(a) (thus the size distribution P(cr)) is permitted to fluctuate, therefore, they can sim- 
ulate a true polydisperse system in the average sense. Comparing with the grand canonical 
ensemble, the semigrand ensemble is more suitable because the total number of particles 
is fixed and the insertion or deletion of particles is not needed. The semigrand ensemble 
is especially more suited to simulate the dense fluid and crystaUi^ 1 ^, which provides us a 
perfect framework to investigate the stability of polydisperse hard sphere crystals. 

In the semigrand ensemble, the particle size distribution P(cr) is not chosen a priori, which 
is obtained only after the simulation has been performecUi 1 ^ 1 ^. This is because the imposed 
physical variables in the simulation are not the composition distribution but the chemical 
potential deference function A/j,(a) = h(<j) — fj,(a r ) (here a r is the diameter of an arbitrarily 
chosen reference component), which is a functional of the composition distribution P(o~). 
Consequently, in order to simulate a system with a prescribed distribution the inverse prob- 
lem A/i(cr) = A / u({P(cr)}) has to be solved. Recently, Escobedo 2 ^ and Wilding et al 21 have 
separately shown that the inverse problem can be solved by a histogram reweighting method. 
Alternatively, a more robust and convenient scheme, the so called nonequilibrium potential 
refinement algorithm(NEPR), was proposed by Wilding^ and works excellent in the grand 
canonical simulation. We will extend the algorithm to the semigrand ensemble, and use the 
extended method to determine the chemical potential deference function A/i({P(a)}) of an 
arbitrarily prescribed composition distribution P(cr) in a semigrand ensemble. The resulting 
forms of Ayu({P(<j)}) are then used to study the stability of the polydisperse hard sphere 
crystals. 

To avoid any confusion with the Helmholtz free energy, we will refer to the free energy 
of the semigrand ensemble as the semigrand free energy. In the semigrand ensemble the 
most stable phase has the lowest semigrand free energy. There are basically two routes to 
follow in the evaluation of the semigrand free energy. One is the thermodynamic integration 
route 2 ^ 1 ^ which determines the free energy of a system by integrating its derivatives along a 
parameter space path connecting the system of interest to a reference system( e.g., Einstein 
solid or ideal gas). The other is the lattice switch method proposed by Bruce et al^, from 
which the free energy difference between monodisperse fee and hep hard-sphere crystals can 
be calculated more directly than the thermodynamic integration method. The method is 



utilized in the canonical and isobaric-isothermal ensemble^. Therefore, We will extend the 
lattice switch approach to the semigrand ensemble in the present work and use it for the 
study of thermodynamic stability of the polydisperse hard sphere system. 

The contents of the remain of this paper are as follows. In section [TT] we formulate the 
statistical mechanics for a polydisperse system within the semigrand ensemble. In section UlTl 
the methodology employed in the work is described and their validity is checked. The com- 
putational details and results are presented in section HVl Finally, we present our conclusions 
in section IVl 

II. THE SEMIGRAND CANONICAL ENSEMBLE 

The most convenient ensemble in the simulation of a polydisperse system is the so called 
semigrand canonical ensembleSCE though other ensembles can also be used. In the SCE the 
total number of particles N and the volume V are fixed while the sizes of each particles can 
be changed. The average particle size distribution is determined by the chemical potential 
difference function A/z(cr). First, lets consider a system of N hard-spheres in a volume 
V, and the distribution of the diameter of the spheres is P{&)- Here we assume that the 
number of particles N is large enough so that the distribution P(cr) can be well defined. 
The Helmholtz free energy of the system is 

A = -PV + N J n(a)P(a)da. (1) 

The semigrand canonical free energy (SCFE) is defined through a Legendre transform^ 

Y = A - N f(fi(a) - fx(a r ))P(a)da. (2) 

Here fi(a r ) is the chemical potential of the reference particle (with diameter o>). SCFE Y 
is a function of temperature T, volume V, total number of particles N and a functional of 
A/x(cr). The partition function for SCE is 
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Here a, and A(cTj) = h/{2^mikT) 1 / 2 are the diameter and the thermal wave length of the 
ith particle, respectively. And Zn is the canonical configuration integral 
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By setting /i ex (o"j) = ^((7j) — /cTln ( — ^ii- J as the excess chemical potential from idea gas, 
the semigrand canonical partition function can be written in a more symmetrical form 

Iff N N 

The SCFE of the system is related to the partition function by the following relation 

Y = -kT\nr(N,V,T,An(a)). (6) 

Thus the stable state can be obtained in the semicanonical ensemble by the minimization of 
the semicanonical free energy, which is the criteria for the stability of the polydisperse hard 
sphere crystal. 

In practical simulation calculations, the diameter of particles is discretised and the cor- 
responding semigrand canonical partition function is 

1 °X ox N 
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where a s and o^are the maximum and minimum values of the particle diameters, respec- 
tively. The above discussion can be extended straight forwardly to the polydispersity of other 
properties of the particles, such as charge dispersity, shape dispersity and mass dispersity 
etc. 



III. THE METHODS 
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Under the semigrand canonical ensemble (SCE), the excess chemical potentials for 
different particles (different diameters) relative to the reference particle(diameter a r ), 
\iex{?i) ~ VexiPr), are given. However, in experiments the fixed quantity is the distribu- 
tion of particle diameters -P(cr), in order to simulate the experimental controllable system 
with given particle size distribution, a proper excess chemical potential has to be chosen 
which can reproduce the required particle size distribution. This is in fact the solution of 
the functional equation Afi ex (a) = A/j, ex ({P(a)}). Wilding has proposed an effective and ro- 
bust procedure, the nonequilibrium potential refinement (NEPR) algorithm 22 , to tackle this 
problem. The algorithm can be used to solve a wide range of the so-called inverse problems 26 
such as obtaining the inter particle interactions from experiment measured structure factors. 
It can also be used in our problem to find the excess chemical potential from the particle size 
distribution. The original NEPR algorithm was developed in the framework of the grand 
canonical ensemble, we extended it to the case of the SCE and used it in the calculation 
reported here. The following is the detailed description of the extension. 

Consider a polydisperse system of hard spheres, the range of the particle diameter is 
cr s — <j\ • • • <Ji ■ ■ ■ a c — (Jl, and the diameter of particles can take c discrete values, <7j, i — 1, 
2, • ■ ■ , c. When c is large enough, the diameter of the particles tends to a continuous variable 
which can resemble the real polydisperse system. The diameter distribution is -P(cr), which 
is normalized in the following way 

c 

5>(*) = 1- (8) 

The excess chemical potential A/j, ex (a) is solved by simulation in a recurrence way. First, 
initial guess of the excess chemical potential is assigned, then it is modified at every few 
Monte Carlo steps according to the instant diameter distribution, P ins (a), which records the 
distribution of particles at the instant of the simulation, the simulation is terminated when 
the average of Pi ns (<j) is the same as the required distribution P(a) within some tolerance. 
Then the A// ea .(cr) is the solution of the problem. The initial value of the excess chemical 



potential is not a vital factor in the calculation process and may be assigned any reasonable 
values, for example, Afi ex (a) = 1 for all diameters. The detail implement is the following. 

1. The particle move 

There are two kinds of particle moves in the simulation, the first kind is the random 
displacement of a randomly chosen particle, which is rejected or accepted depends on whether 
the new position overlaps to other particles or not, the second kind is the expansion or 
retraction of a randomly chosen particle, which is named as breathing move in literature^, 
the probability of acceptance of a breathing move which does not result in an overlap is 

P acc = min |l,exp{/5(A/i ex (a-) - Afi ex (ai))} > , (9) 

where <Tj and a i are diameters of the zth particle before and after the test move. The move 
is rejected if it results in an overlap with other particles. 

2. The iteration 

For a given particle size distribution, the excess chemical potential is calculated by a 
Monte Carlo iteration procedure. The central quantity in this procedure is the instantaneous 
particle size distribution P ins (a), which is the histogram of the particle size distribution at the 
instant of the simulation and updated during the simulation. Another important quantity 
is the average particle size distribution P(cr), which is the average of the instantaneous 
particle size distribution in the simulation. The excess chemical potential is updated by the 
Wilding's scheme^ for every short intervals. The Wilding's scheme in this iteration is given 
by 

A^(a) = AaU<0 ~ H ( Pinsi p 7 f {<T) ) Va. (10) 

Here P{<j) is the given particle size distribution, 7; is a modification factor of the ith iteration. 
For a given modification factor, the average size distribution P(o~) is also recorded during 
the simulation. When the difference of the average size distribution and the given particel 



size distribution is less than a specified value £ 
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one loop of the iteration is finished. The modification factor is then reduced by a factor 
1/n where n is a small integer, j i+ i = ji/n, and the excess chemical potential of the last 
iteration is used as the initial input and start the next iteration. The iteration continues 
till the modification factor 7 reaches a very small value, and the resulted excess chemical 
potential is then regarded as the solution of the problem. In practical calculations, the 
convergent criteria for 7 , the threshold £ and the reduced factor n are tuned to reach both 
high efficiency and accuracy. 

The NEPR algorithm for SCE was tested by the simulation of a polydisperse hard sphere 
fluid and a polydisperse hard sphere crystal with fee structure. In the tests the particle 
size distribution is chosen to be the Schultz distribution, which is the most studied model 
distribution in polydisperse systems. The distribution is 



P(a) 



1 f z + 1 
z\ V a 



2+1 



a exp 



z + 1 



a 



a 



(12) 



where a is the average diameter of the particles and z controls the width of the distribution. 
In the Schultz distribution, the range of the diameter is [0, +00), however, in a simulation 
calculation with finite number of particles, cuts off of the up and the lower limit of the 
distribution may be specified for convenience of computation. In the test studies the effect 
of the cutoff is not studied, the emphases is on the effect of polydispersity to the physical 
properties of the system. On the other hand, small particles may enter into the interstitial 
space of crystals and induce instabilities of crystal structure, this is beyond the subject of 
this study though it is an interesting subject of research. 

In the test simulation, the threshold £ = 0.15, the initial modification factor 70 = 0.01 and 
the reduce factor is 1/2. The termination criteria is 7 < 0.0001. For the hard sphere liquid, 
the volume fraction <fi = 0.3 and the dispersity 5 = 14.2%; for the fee hard sphere crystal, 
the volume fraction <fi = 0.6 and the dispersity 5 = 3.8%. Figure Q] and figure [2] are the 
simulation results. Figure QJa) is the calculated excess chemical potential for the truncated 
Schultz distribution as function of the diameter of particles in the liquid state, figure Wfo) is 



the comparison of the given truncated Schultz distribution and the distribution generated 
with calculated excess chemical potential, the agreement between the two is excellent. Figure 
|2(a) and figure [2(b) are the calculated excess chemical potential and the comparison of given 
and generated distributions in the crystal state, respectively, the agreement is also excellent 
as in the liquid case. 

B. Lattice switch 

The free energy difference between the fee and the hep hard sphere crystals is extremely 
small, in order to obtain a reliable result for the difference, we need to find an accuracy 
method of calculation. There are different methods suggested in the past in the studies of 
monodisperse hard sphere crystals and many results were obtained for the problem. The 
lattice switch method(LW)- 1 ^ developed recently is a high precession method in the calcu- 
lation of the free energy difference. The method has extended successfully to the calculation 
of the liquid-solid transition of monodisperse hard sphere systems^ and also used in the 
studies of soft sphere systems^*^. In this subsection, we extend it to the SCE. 

A detailed presentation of the LW method for monodisperse hard sphere crystals in the 
canonical ensemble can be found in references' 1 ^. Here we give a quick sketch of the method 
in the context of polydisperse hard sphere crystals. The system contains N hard spheres in 
volume V with periodic boundary conditions. The hard spheres can be in the fee or the hep 
structures respectively, and the spatial positions of the particles are specified by the position 
vectors Yf cci or Th cpi for the ith particle in the fee structure or hep structure, respectively. 
The position vectors can be decomposed as 

where the subscript a may represent fee or hep. The R/ CC j and Hhcpi are lattice vectors of 
the idea structures, Uf cci and u hcpi are the displacements from the idea structures. 

In principle, the displacement can be any vectors that only constrained by the geome- 
try of the simulation box, however, in the crystal phase with dispersity smaller than the 
terminal dispersity, the displacements are naturally cutoff in the simulation time scale. We 
use {u, a} to represent the phase space of the polydisperse system, here a denotes the di- 



ameters of N particles. Each structure a(a = fee or hep) associates a set of displacements 
{u, cr} a - In a typical simulation, in which a representative subset of the displacements for 
one structure is sampled, the transition from one structure to another can not happen be- 
cause the transition probability between structures is extremely small. The spirit of the 
lattice switch method is that switch the ideal lattice vectors from one structure to another 
while keep the displacements frozen. The two sets {u, a}/ cc and {u, cr}hcp have a common 
intersection {u,a}f cc f]{u,a}h cp , which provides a gate to relate the two structures. All 
allowed(non overlap) configurations accessible by simulation which are associated with fee 
and hep structures can be divided into three subsets, (a), all the displacements allowed by 
the fee structure but not allowed by the hep structure, which we denote as {u, o-}f cc _ hcp , 
(b). all the displacements allowed by the hep structure but not allowed by the fee structure, 
which we denote as {u,a}h cp -fcc, (c). the displacements allowed by both the fee structure 
and the hep structure, denoted as {u, o-}f cc (^ hcp . 

The semigrand canonical ensemble partition function of the two structures can be written 
as 

1 „ AT „ N N 

r («) = An \3N( \ / Y[ dcTi / YldUiXexp{(3(^2A^ ex (ai)-(f)(u,a))}, (14) 
iV.A {cr r ) j ae{a} . =i Jue{u} i=1 i=1 

where is the potential energy of the system, which is oo or for the hard sphere system. 
The SCFE for the structure a is 

Y(a) = -kT\nT(a). (15) 

The SCFE difference of the two structures can be written as 

nice) 



Yhcp — Yfcc — kT hi 
= kTlu 



T(hcp) 

Pfcc 
±hcp 

fec- 
-L hep— fee T -Lfcc{~}hcp 



kT dn f cc ~ hcp — "f cc C\ hc P Qg\ 



Where Pf cc _hc P , Phcp-fcc and Pf cc (~\hcp represent the probabilities of three subsets 
{u, cr}f cc -hcp, {u, cr}hc P -fcc and {u, cr} fccrihcp respectively 
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It is clear from the above discussion that the calculation of semigrand free energy dif- 
ference from simulation may be proceed as follows: first, the excess chemical potentials 
for both structures are determined by the iteration procedure described in the last section. 
Then starting from any structure, particle size distribution and displacement, regard all 
of the displacement, particle size and the structure index as random variables, and make 
Monte Carlo moves. In principle, the system will move among the fee — hep, hep — fee 
and fccf]hcp states, by recording the number of microstates corresponding to the three 
macrostates, the probabilities of each macrostate can be obtained and then follows the free 
energy difference. However, this prescription is not practical in real simulations, the system 
will trap in either fee — hep or hep — fee macrostate in the simulation period simply because 
the number of microstates of fee — hep and hep — fee are much larger than the number of 
microstates of fee f] hep. In order to overcome this difficulty, the bias sampling can be used. 
To achieve this goal, we first define an order parameter ^(u, a) for the displacement field, 

^#(u, a) = M(u, a, hep) - M(u, a, fee). (17) 

Here M(u, a, hep) and M(u, a, fee) represent the number of overlap pairs of the hep and 
fee structure for all samples of the displacements. The order parameter ^ can take values 
of 0, ±1, ±2, • • • , where ^ = corresponding to the macrostate feef] hep. In the simu- 
lation process one of the M(u, a, hep) and M(u, a, fee) has to be zero since the domain of 
random walk is {u, a}/ cc |J{u, a}hc P - The free energy difference can be represented by the 
macroscopic order parameter <M as 



Y hcp -Y fcc = kT ^^- . (18) 



here P{^) is the probability that the order parameter takes value ^#. Now we regard 
each value of the order parameter corresponding to a macroscopic state, biased sample 
alorithm^^i^ is to sample the system according to a probability so that the rate of visits 
to every macrostate is basically the same. If the sampling probability is the inverse of P(^), 
then the rate of visits to each macrostate is exactly the same. Unfortunately, P(^) is the 
quantity we are looking for which is unknown before calculation. The problem was solved 
by several different methods in the context of density of states calculations, from which the 
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multicanonical method 31 and Wang-Landau method 3 ^ are the most used. These methods, 
when used in the current problem, amount to starting the calculation with an initial guess 
of the probability P(^#), sampling the system with inverse of this probability, and modify 
P(^C) according to the visit rates to the macrostates till the visit rates are constant for 
all macrostates within a given tolerance, then the resulted P(*4f) and visit rates together 
will give an accurate estimate of the free energy difference. One of the implementaion is to 
introduce a weight function r\{^M\ sampling the system with e^^* Afj.ex{o-i)-<f>(u))+ri(jf) ^ ca j cu _ 
lating the probability distribution of the order parameter ^#, P r? (^#), modifying the weight 
function r]{^) till P n {^£) is close to constant and then obtain the required probability 
P(e/#) through 

P(Jt) = P v {J?)e~^\ (19) 

Based on these discussions, the acceptance probability of the sampling is summarized in the 
following expression 



j i f , 



PacA-M, u, a -»• Ji , u , a = mm <M , — — — — — -— \ , 20 

where(^#, u, a) and {^M , u , a ) are the order parameters, displacement fields and diameters 
of the system corresponding to states before and after the test move, respectively. 

IV. COMPUTATIONAL DETAILS AND RESULTS 

In this section we use the extended NEPR and lattice switch method described in the last 
section to study the stability of the polydisperse hard sphere crystal. The distribution of the 
hard spheres is chosen to be the Schultz distribution as was used by many researchers. The 
chemical potential difference obtained from this distribution with fee lattice can reproduce 
fairly accurate Schultz distribution for the specified hep lattice, this is because of the differ- 
ence between the two structures is very small(see Fig. [3]). In fact, if we specify a chemical 
potential difference, we may produce the same particle size distributions in both of the fee 
and the hep structures within the statistical errors. It is well known that there is a terminal 
poly-dispersity in polydisperse crystals above which the bulk crystal may not stably existed. 
The terminal polydispersity St obtained from recent simulation is about 5 ~ Q%r^^. In 
our simulation we set the maximum dispersity to be 4% so that the stable crystal can be 
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simulated. The simulation boxes are set up to suit the idea fee and hep crystal, periodic 
boundary conditions are used in the calculation, the initial configuration are idea fee and 
hep lattices, respectively. The Wang-Landau sampling method was used to obtain a crude 
estimation of the weight function of the order parameter, then the multicanonical algorithm 
is used to refine the result. The calculated results are shown in Table HJ for the polydis- 
persity used in the simulation, the fee lattice has the smaller free energy and more stable 
than the hep crystal. In the case of 5 = 4%, there are possibilities that the smallest sphere 
may jump from the cage of its lattice positions so that a defect may be created, to avoid 
this situation we used a larger volume fraction as shown in the last column of Table [H The 
effect of finite size effect was studied in the case of medium polydispersity by enlarging the 
simulation boxes with fixed volume fraction, it was found that the value of the free energy 
is affected by the box size slightly but the relative stability is unchanged. 

In order to have a clear picture of the relative stability of the structures, we plotted the 
probability distribution of the macrostates for two different polydispersities in figures 0H5l 
For convenience of comparison, we plotted the distribution of fee and hep in the same half 
plane, in the hep case the absolute value of the order parameter is used. From the figures we 
see that there is a maximum of probability for each structure, and the probability maximum 
of fee structure is larger then that of the hep structure which means that fee structure is 
more favorable. Figure [6] shows that the ratio of the probability of "gate" states [ytf = 0) 
to the probability maximum is about 1CT 35 which means that the "gate" states will never 
be reached if a simple sampling scheme is used. Considering the extreme difference of the 
probability between the "gate" state and the maximum state, the refinement of macrostates 
and bias sampling are the necessary scheme to obtain meaningful results. 

The particle size distribution used in the calculation is the truncated schultz distribution, 
this is the widely used distribution in theoretical studies. The results based on this distribu- 
tion may not be extended to all polydisperse systems, however, we expect that it represents 
a class of polydisperse systems. The distribution dependence of structures of polydisperse 
systems requires detailed computations of various polydisperse systems. It should also be 
noted that the accuracy of our calculation is limited by the computation resources, the num- 
ber of particles is pretty small, typical Monte Carlo steps are 50 million which is still too 
small to determine accurately the relationship between the difference of semigrand canonical 
free energy and the polydispersity. On the other hand, the conclusion that the fee is more 
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favorable than the hep structure for the polydisperse hard sphere crystals with truncated 
Schultz distribution of diameters of spheres is clear and reliable. 

V. CONCLUSION AND DISCUSSION 

Polydispersity of colloid system is common in artificial colloids. We studied here the effect 
of polydispersity on the stability of structures of colloid crystals, and found that fee structure 
is more stable than the hep structure, which is the same as the monodisperse case with the 
same calculations. To study the problem we have extended the NEPR algorithm and the 
lattice switch method to the semigrand ensemble. The extension provides a powerful tool in 
the studies of other thermodynamical problems of polydisperse systems. A direct application 
is the determination of the phase diagrams of the polydisperse systems which may replace or 
at least complement the current Gibbs-Duhem integration method^ 3 ^. The monodisperse 
system has already studied by the original LW method 25 . The above extension can also be 
extended to the problems of soft sphere system simulations by some extra techniques' 3 ^. 

Acknowledgments 

The work is supported by the National Natural Science Foundation of China under grant 
No. 10334020 and in part by the National Minister of Education Program for Changjiang 
Scholars and Innovative Research Team in University. 



Electronic address: hrma@sjtu.edu.cn 



1 B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 (1957) 

2 W.W. Wood and J. D. Jacobson, J. Chem. Phys. 27, 1207 (1957). 

3 W.G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 (1968). 

4 B. J. Alder, B. P. Carter, and D. A. Young, Phys. Rev. 183, 831 (1969). 

5 L. V. Woodcock, Nature (London) 384, 141 (1997). 

6 P.G. Bolhuis, D. Frenkel, Siun-Chuon Mau and David A. Huse, Nature 388, 235(1997). 

7 A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev. Lett. 79, 3002 (1997). 

8 S. Pronk and D. Frenkel, Phys. Rev. Lett. 90, 255501 (2003). 

14 



9 S. K. Kwak and D. A. Kofke, J. Chem. Phys. 122, 176101 (2005) 

10 P. N. Pusey and W. van Megen, Nature (London) 320, 340 (1986). 

11 PG. Bolhuis and D. A. Kofke, Phys. Rev. E 54, 634 (1996). 

12 M. Fasolo and P. Sollich, Phys. Rev. Lett. 91, 068301 (2003). 

13 S. Phan, W. B. Russel, Z. Cheng, J. Zhu, P. M. Chaikin, J. H. Dunsmuir, and R. H. Ottewill, 
Phys. Rev. E 54, 6633 (1996). 

14 S. Phan, W. B. Russel, J. Zhu and P. M. Chaikin, J. Chem. Phys. 108, 9789 (1998). 

15 S. Martin, G. Bryant, and W. van Megen, Phys. Rev. E 67, 061405 (2003). 

16 H. J. Schope, G. Bryant and W. van Megen, Phys. Rev. Lett. 96, 175701 (2006). 

17 J. G. Briano and E. D. Glandt, J. Chem. Phys. 80, 3336 (1984). 

18 D. A. Kofke and E. D. Glandt, J. Chem. Phys. 87, 4881 (1987). 

19 M. A. Bates and D. Frenkel, J. Chem. Phys. 109,6193 (1998). 

20 F. Escobedo, J. Chem. Phys. 115, 5642 (2001); 115, 5653 (2001). 

21 N. B. Wilding and P. Sollich, J. Chem. Phys. 116, 7116 (2002). 

22 N. B. Wilding, J. Chem. Phys. 119, 12163 (2003). 

23 D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 (1984) 

24 D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, San Diego, 1996). 

25 N. B. Wilding and A. D. Bruce, Phys. Rev. Lett. 85, 5138 (2000). 

26 D. Levesque, J. J. Weis, and L. Reatto, Phys. Rev. Lett. 54, 451 (1985). 

27 M. R. Stapleton and D. J. Tildesley, J. Chem. Phys. 92, 4456 (1990) 

28 A. D. Bruce, A. N. Jackson, G. J. Ackland, and N. B. Wilding, Phys. Rev. E 61, 906 (2000) 

29 A. N. Jackson, A. D. Bruce, and G. J. Ackland, Phys. Rev. E 65, 036710 (2002) 

30 J. R. Errington, J. Chem. Phys. 120,3130 (2004) 

31 B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992) 

32 F. Wang, D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001); Phys. Rev. E 64, 056101 (2001). 

33 J. S. Wang and R. H. Swendsen, J. Stat. Phys. 106, 245 (2002) 

34 D. A. Kofke and P.G. Bolhuis, Phys. Rev. E 59, 618 (1999) 



15 






0.576 


0.576 


0.576 


0.576 


0.576 


0.602 


5 








1.1% 


2.5% 


2.5% 


4% 


N 


216 


1728 


216 


216 


1728 


216 


A/ x 10 5 


133(3) 


113(3) 


139(11) 


133(16) 


110(15) 


170(20) 



TABLE I: Parameters and simulation results: (j) is the volume fraction of the system, 5 is the 
polydispersity of particle diameters (same for both fee and hep structures), N is the total number 
of particles of the system, Af = (Yh cp — Yf cc )/N is the semigrand canonical free energy difference 
between hep and fee structures, in unit of kT. The numbers in the parentheses is the uncertainty 
of the result. The results of monodisperse hard sphere system(from reference^) are listed in the 
first two columns for comparison. 
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FIG. l:(a), The solved excess chemical potential difference of hard spheres as function 
of particle diameter a in the fluid state, (b), The line is the plot of the Schultz function, 
and the dots are the particle diameter distribution obtained from simulation by using the 
Afi ex (a) plotted in (a). 

FIG. 2: The same as figure [H for fee solids. 

FIG. 3: The calculated excess chemical potentials for the truncated Schultz distribution 
of particle diameters from NEPR method. The line is for the fee structure and the points 
are for the hep structure, the difference is smaller then the statistical errors. 

FIG. 4: Left: the excess chemical potential as function of the particle diameter, insert 
is the particle size distribution for fee crystal. Right:the distribution of macrostates (order 
parameters) for fee (open circles) and hep (filled square) crystals. The polydispersity 
5 = 1.1%, volume fraction (ft = 0.576. 

FIG. 5: Same as figure IH 5 = 4%, volume fraction <p = 0.602. 

FIG. 6: The logarithm of the probability distribution of macrostates of hard sphere 
crystals, the solid line is for the fee crystal and the dashed line is for the hep crystal. 
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